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I .  INTRODUCTION 


The  K-integrator  for  stiff  ordinary  differential  equations  was 
developed  in  the  late  1960 's  by  one  of  us  (MDK) ,  and  has  been  repeatedly 
improved  since  that  time^"^.  It  is  intended  to  solve  systems  of  the  form 


Y(t)  =  G  (Y (t) ,  t) ,  te  [a,b] ,  (1) 

YCa)  =  Y0  , 

where  Y(t)  =  dY(t)/dt, 

Y(t)  =  [Y1 (t) ,  Y2 (t) ,  ...  ,  YN(t)]T  , 
and  G (u, t)  =  [G1 (u,t) ,  G2(u,t),  ...  ,  GN(u,t)]T. 


The  concept  of  stiffness  is  difficult  to  define  formally.  However, 
it  can  be  described  in  terms  of  the  Jacobian  matrix  J(u,t),  where  the 
element  in  row  i  and  column  j  is 

Jlj(u,t)  =  9G1(u,t)/9y-*  .  (2) 


A  stiff  system  will  have  one  or  more  eigenvalues  in  the  Jacobian  whose 
real  parts  are  negative  and  large  in  modulus.  As  a  result  the  corre¬ 
sponding  components  in  the  solution  will  decay  very  rapidly  in  comparison 
to  the  other  terms  present. 


M .  D .  Kregel  and  E .  L .  Lortie^  " Description  and  Comparison  of  the 
K-Method  for  Performing  Numerical  Integration  of  Stiff  Ordinary 
Differential  Equations" 3  BEL  Report  No.  1722,  July  1974  (ADA  #A 003855). 


2 

M.  D.  Kregel  and  J.  M.  Heimerl,  "Comments  on  the  Solutions  of  Coupled 
Stiff  Differential  Equations",  BEL  Memorandum  Report  No.  2769,  July  1977; 
or  Proceedings  of  the  1977  Army  Numerical  Analysis  and  Computers 
Conference j  AEO  Report  77-2,  November  1977 ,  pp.  553-562  (ADA  #A043122) . 

$ 

T.  P.  Coffee ,  J.  M.  Heimerl,  and  M.  D.  Kregel,  "A  Numerical  Method  for 
Large  Stiff  Systems  of  Ordinary  Equations" ,  Transactions  of  the  24th 
Conference  of  Army  Mathematicians,  ARO  Report  79-1,  January  1979. 
pp.  249-257. 
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It  will  be  useful  to  rewrite  the  system  (1)  using  the  diagonal 
Jacobian  elements.  We  define  the  diagonal  matrix  R(Y(t),t)  and  the 
vector  F(Y(t),t)  with  components, 


R1 (Y  (t) ,t)  =  -J11 (Y (t) , t) ,  and  (3) 

F1  (Y  (t)  ,  t)  =  Y1  (t)  +  R1(Y(t),t)Y1(t). 


Then  equation  (1)  becomes 

Y  (t)  =  F(Y(t),t)  -  R(Y  (t)  ,  t)Y(t)  .  (4) 

The  term  R1  measures  the  sensitivity  of  Y1  with  respect  to  changes  in 

Yi . 


The  system  is  solved  at  a  set  of  discrete  points  in  time.  Thus, 

n 


step  sizes,  h. ,  and  times,  t  =  a+  £  h.,  are  introduced.  Approximations 

i=l 

y^  to  Y (t^)  are  then  produced  for  n=l,2,  ...  by  the  integrator. 

These  approximations  will  be  found  using  multistep  formulas  of  the 


form  y 
yn 


q 


i 

i=l 


a . 
1 


y 


n-i 


Bi 


n-i 


(5) 


The  truncation  error  will  depend  on  the  order  of  the  formula.  For  stiff 
equations,  there  is  also  the  problem  of  stability.  That  is,  a  small 
error  in  one  step  may  propagate  and  grow  in  subsequent  steps .  Standard 
non-stiff  methods  will  be  restricted  to  very  small  step  sizes  to  preserve 
stability.  Dahlquist^  has  shown  that  to  maintain  reasonable  stability 
for  stiff  systems,  the  formulas  must  be  implicit,  that  is,  3^  f-  0. 

Most  stiff  multi-step  methods  use  a  modified  Newton  iteration  to 
solve  (5).  First  an  explicit  predictor  yn^  is  found.  Then  its  time 
derivative  and  corresponding  corrector  are  defined: 


,t  ) 
5  w 


,t  )y 
n  'n 


and 


(6) 


^ G .  G .  Dahlquist ,  AMS  Symp.  Appl.  Math .  f 1963 ) ^  147 . 
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+ 


(7) 


q 


l 

i=l 


a .  y 
l'n-i 


q 

l 

i=l 


g.y  n  +  h  g  y 
i'n-1  n  crn 


Finally  the  vector  A 

n 


is  defined  such  that 


y 


n 


y  +  A  . 
n  n 


(8) 


P  c 

An  approximation  dn  to  An  can  be  found  from  the  difference  (yn  -yn  ) , 
by  expanding  yn  in  a  series  about  the  ynp  and  truncating.  This  procedure 
leaves  the  following  system  of  linear  equations  to  be  solved. 


[h3  J 
1  o  n 


(9) 


p 

where  Jn  =  J (yn  ,tn)  and  I  is  the  identity  matrix.  The  accuracy  of  the 
approximation  to  yn  can  be  found  by  monitoring  the  size  of  the  quantity 
(y  P2  -  y  C2)  where 
wn  'n  J> 


P2 


y 


n 


+  d 

n 


and 


(10) 


q 

l 

i=l 


a.  y  .  + 
iy  n-i 


n 


q 

l 

i=l 


$.y 

l'n- 


h  3  y 
n  cr  n 


P2 


(11) 


If  necessary,  the  Newton  iteration  can  be  repeated. 

3 

Solving  the  system  of  equations  (9)  requires  approximately  N  /3 
multiplications  and  divisions,  where  N  is  the  number  of  equations.  Thus, 
in  general  the  larger  N,  the  greater  the  computation  time.  Reduction  of 
this  time  has  been  an  objective  of  several  algorithms  including  the  one 
presented  herein. 

One  approach  to  this  problem  is  the  well-known  algorithm  DIFSUB,  by 
C.  W.  Gear^.  This  is  a  variable-order  method,  based  on  formulas  of  the 
form 


5C.  W .  Gear ^  " Numerical  Initial  Value  Problems  in  Ordinary  differential 
Equations” >  Prentice-Hall ,  1971 . 
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n 


q 

=  l 

i=l 


a .  y 

1  n-i 


+  h  60  V 


1  <  q  <  6. 


(12) 


The  program  automatically  uses  the  higher  order  formulas  if  more  accuracy 
is  required.  All  the  formulas  are  based  on  a  fixed-step  size  h.  To 
change  the  step  size,  the  appropriate  values  for  y  are  found  by 
interpolation. 

For  a  given  step,  the  matrix  based  on  the  Jacobian  is  found  and 
inverted.  The  program  reevaluates  this  matrix  only  when  it  fails  to 
obtain  convergence  in  three  Newton  iterations.  One  therefore  does  not 
have  to  invert  a  matrix  every  step. 

More  recently,  variable-step  formulas  have  been  developed.  Any 
sequence  of  step  sizes  h^  can  be  used.  The  coefficients  and  3-j_  in 
(5)  are  variable  and  depend  on  the  previous  step  sizes.  In  practice, 
such  codes  have  been  more  stable. 

One  such  code  is  EPISODE,  developed  by  Byrne  and  Hindmarsh^.  It 
uses  a  variable-order,  variable-step  formula.  Like  DIFSUB,  it  uses  an 
aged  Jacobian.  An  LU  decomposition  is  computed  instead  of  a  matrix 
inversion. 


The  K-integrator  uses  a  fixed-order,  variable-step  formula.  Except 
for  starting,  it  always  uses  a  third-order  formula.  The  lack  of  higher 
order  formulas  will  not  be  important,  provided  only  moderate  accuracy 
is  required. 


The  major  innovation  in  the  K-integrator  is  the  method  of  solving 
the  set  of  equations  (9).  Both  DIFSUB  and  EPISODE  approximate  (h30Jn-I) 
by  the  values  of  the  matrix  at  an  earlier  time  step  tn_-[.  The  K-integrator 
evaluates  the  matrix  at  each  time  step;  but,  if  possible,  it  is  one  of 
reduced  size.  Finding  a  suitable  approximation  to  any  of  the  dn^ 
permits  the  corresponding  row  and  column  to  be  eliminated  from  the 
matrix. 


To  illustrate,  we  write  (9)  in  component  form; 


i.e.  , 


Gk 


=  hfi 


N 

l 

j=i 


kj 


n 


[he  J 
L  o  n 


kk 


-  l]d 


(13) 


D.  Byrne  and  A.  C. 


Hindmarsh,  ACM  Trans_.  Math.  Software  (1975),  71. 


O 

O 


We  attempt  to  approximate  d^K  by 


dnk  -  (XnPk  -  ^nCk^hVnkk  "  « 


(14) 


There  are  several  cases  where  (14)  is  a  reasonable  approximation. 
First,  if  h  is  very  small,  such  as  at  the  beginning  of  the  integration, 

so  that  |h3Q  £  ^ |  <<  |d  k|,  then  the  off-diagonal  terms  may  be 

neglected,  Second,  even  if  h  is  large,  the  equation  for  yn  may  be 

weakly  coupled  to  the  other  equations,  i.e.,  |j  kkd  k|  >>  I  Y  J  ^  d  ^ I 

n  n  1  1 .H,  n  n  1 

•  Pk  rv 

Finally,  at  any  value  of  h,  yn  may  be  a  good  approximation  to  y^  , 

and  so  |y  ^  -  y  ~  0  ~  d  ^ . 

1  n  n  1  n 

In  general,  the  set  of  equations  (13)  is  solved  in  two  stages. 

First  a  set  of  indices,  L,  is  determined  for  which  (14)  is  valid,  i.e.. 


a  1  r  p£  CJL  /ru0  ,  M 

dn  ■  %  ~  yn  >/tWo  Jn  '  l>  >  UL' 


(15) 


Then,  using  these  values,  the  system  (13)  becomes 


y  Pk  -  y  Ck  -  he  l  J  KJ  d  1  =  hd  !  J  Kld  1  .  [k»  J  lid  \ 
n  Jn  o  >T  n  n  o  .5.  n  n  L  o  n  J  n  * 

jeL  jtL 

k  i  L.  (16) 


kJ  a  j  - 


kj  A  J 


kk 


The  reduced  system  (16)  is  then  solved  using  an  LU  decomposition  and 
back  substitution.*  Where  the  reduced  system  involves  a  small  number 
of  the  original  set  of  equations,  for  example  when  h  is  small,  there 
will  be  a  significant  reduction  in  computation  time. 

The  difficulty  in  using  this  algorithm  lies  in  determining  which 
set  of  equations  can  be  solved  with  sufficient  accuracy  using  (14) . 

The  version  reported  here  monitors  the  error  terms  |yn^ ^  _  ynC2k|  anj 
uses  them  to  determine  the  set  of  indices  L  for  the  next  time  step. 

The  details  are  given  in  the  next  section. 


* Earlier  versions  of  the  K-method  (Refs.  1  and  2)  included  a  Gauss- 
Seidel  iteration  procedure  before  the  matrix  factorization.  The  version 
reported  here  omits  this  iteration  procedure. 
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II.  THE  NUMERICAL  METHOD 


The  K-integrator  is  based  on  a  fixed  third  order  formula  of  the 

form 


ai  y  i  + 

1  yn-l 


'2  yn-2  *  Vn-3+h  (60  >n 


+  61  Vl  * 


®2  V2> 


It  is  a  variable  step  formula.  If  h,  h',  and  h"  are  the  last  three  step 
sizes,  then 


e2  =  0.05 
=  0.35 
p  =  h'/h 

q  =  Ch '  +  h")  /h 

80  =  [(1+P)  (1+q)  -  &2P(p-q)  "  6jP<l]/C3+2p+2q+pq) 

a3  =  [ l+p-e2P2  "  e0 C3+2p) ] /  [q2 (p-q) ] 

a2  =  (B2  +  -  a3q  +  -  l)/p 

al  =  1  "  a2  "  a3  * 

For  constant  step  size,  the  stability  can  be  described  in  terms  of 

parameters  a  and  D.  A  formula  is  said  to  be  A  (a,  D)  stable,  ae  (0,  2.), 

if  all  numerical  solutions  to  Y  =  -AY  converge  to  zero  as  n  -*■  °°  with  h 
fixed  for  all  |  arg  (-  Ah)  |  <  a,  |  A  |  i  0  and  for  all  Re  (hA)  <_  D7. 

This  combines  the  features  of  the  A (a)  -  stability  of  Widlund8  and  the 
stiff  stability  of  Gear8.  For  the  above  formula  a  =  81°  and  D  =  0.4. 

The  truncation  error,  neglecting  higher  order  terms,  is  -  0.0966  h^Y(4)  (t) . 
If  the  step  size  is  not  held  constant,  the  truncation  error  will  vary 
slightly. 


?G.  K.  Gupta ,  Math,  of  Comp.  30  (1976),  417. 

80.  B.  Widlund,  BIT  7  (1967),  65. 
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In  practice  the  K-integrator  uses  modified  forms  of  equations  (9) 
and  (10)  that  are  somewhat  more  convenient  for  chemical  kinetics  problems, 
where  it  has  mainly  been  applied.  The  modified  form  of  equation  (9)  is 


A  g 
nsn 


(18) 


where  A  ^  =  [hg  J  ^  -  6.  .ly  and 

n  L  o  n  ij  n 


1  i  !/  Pi 
g  =  d  /y 
&n  n  J  n 


In  this  notation  equation  (10)  becomes 


P2i  Pi,.,  i^  xt 

yn  =  yn  t1  +  gn  ^  »  i  ■  1,2,...,N.  (19) 


The  rationale  for  (18)  and  (19)  is  the  following.  For  a  given  network 
or  set  of  chemical  reactions  at  a  known  temperature,  the  rate  at  which 
each  reaction  proceeds  equals  a  constant  multiplied  by  the  concentration(s) 
of  the  chemical  species  involved.  The  modified  Jacobian,  with  components 
J^Jy^i,  can  be  generated  by  adding  and  subtracting  rates,  which  is 
computationally  easy  to  do.  (To  recover  the  Jacobian  matrix  written  in 
terms  of  the  rates  alone,  we  divide  by  the  yn11-)  In  the  modified  nota¬ 
tion,  the  diagonal  approximation  (14)  becomes 


Jn 


(y. 


Pk 


n 


-  y. 


Ck 


n 


)/A. 


kk 


n 


(20) 


In  order  to  minimize  the  computation  time,  we  find  as  many  of  the 
gn  as  possible  using  (20).  As  discussed  in  the  Introduction,  it  is 
expedient  to  make  |ynPk  -  yn^l  as  small  as  possible.  To  this  end  we 
have  developed  an  unorthodox  predictor.  While  not  essential  to  the  code, 
it  does  increase  its  efficiency. 

The  usual  predictor  is  obtained  from  an  explicit  multistep  formula. 
But  difficulties  may  arise  when  y£  is  used  in  equation  (8)  to  obtain 
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Since  the  system  is  stiff,  the  quantities  R1(y^,  tn)  from  equation 
(4)  will  often  differ  in  absolute  value  by  many  orders  of  magnitude.  In 
order  to  finish  the  integration  in  a  reasonable  time,  h  must  become 
large  enough  so  that  some  of  the  terms  hPoR^  become  very  large.  This 
means  small  errors  in  y^  can  become  greatly  magnified  in  y^. 

This  effect  can  be  reduced  by  considering  the  form  of  equation 
(17)  for  the  final  value  y  .  Using  (4)  this  can  be  written  as 


y. 


n  =  al  Vl  +  *2  yn-2  +  “5  Vs  +  h  31  Vl 


(21) 


* h  62  yn-2  *  h  60  F^„.  V  -  h  60R<V  V  V 

We  shall  derive  a  formula  for  y  ^  that  follows  this  general  form. 

First,  we  obtain  predicted  values  for  F  and  R,  using  the  formulas 


Fn  “  T1  Fn-1  +  y2  Fn-2  *  ^3  Fn-3 


(22) 


R  =  Yi  R  ,  +  y,  K  0  +  Yt  R 

n  1  n-1  2  n-2  1 3  n- 


3  - 


The  quantities  yj»  Y2>  an<3  Y3  are  determined  uniquely  for  any  step  sizes 
by  imposing  the  condition  that  the  formulas  be  second  order.  For  the 

(3) 

(y,  t)  and 


case  of  constant  step  size,  the  truncation  errors  are  h  F 
(y,t) ,  neglecting  higher  order  terms. 

Substituting  these  predicted  values  into  equation  (21)  ,  we  obtain 


,  (3) 

h  R 


a.  y  -  +  cu  y  „  +  a„y  _ 
L  7n-l  2  /n-2  3  'n-3 


(23) 


+  hg.y  .  +  h  L  y  „  +  h  FP  -  h  g~RP  yP  . 

1  'n-1  2  Jn-2  On  0  n  /n 


This  equation  can  be  solved  for  y  . 
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Neglecting  higher  order  terms. 


y 


i 

n 


4  (3]i  ^3)i 

h4  f$0  (  F^J1(y,  t)  -  R 


(y>  t)  yPi) 


1  +  h  30  Rhyn,  tn) 


(24) 


i  Pi 

As  h  &q  R  becomes  large,  the  accuracy  ofyn  tends  to  improve.  Of 

course,  errors  still  occur  in  and  R^1;  but,  we  no  longer  obtain 

large  errors  just  because  h  Bq  R1  is  large.  In  fact,  in  the  special 

case  where  F  and  R  are  constant,  yn  -  y£  =0. 

After  the  first  few  steps,  a  further  heuristic  modification  is 
made  in  yP  .  We  define  an  error  term  for  the  predictor  by  the  equation 


Cy, 


P2i 


Pi,  ,  Pi 

Yr,  )/y« 


(25) 


where  yn  is  the  final  accepted  predicted  value.  Then  we  define  a 
weighted  average  of  these  errors  by 


w  1  =  0.5  E  Pl  +  0.5  W1  ,,  (26) 

n  n  n-1 


Finally  in  the  next  step  we  define  a  modified  predictor  whose  components 
are  given  by 


?Mi 

yn+l 


Pi 

Vl 


(1  + 


w  )  . 
n 


(27) 


This  modified  predictor  is  used  in  equation  (7)  to  obtain  the  corrector 


The  purpose  of  this  modification  is  to  detect  any  systematic  errors 

in  the  predictor.  We  assume  that  the  accuracy  of  y£  is  similar  to  the 

accuracy  of  yP_^.  If  the  errors  are  not  systematic,  the  fact  that  (26) 

is  a  weighted  average  will  tend  to  minimize  any  errors  introduced  by 

the  w  1 . 
n 
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The  existence  of  systematic  errors  in  an  explicit  predictor  can  be 
shown  for  the  simple  case  of  a  single  equation  Y  =  -  AY,  Y  (a)  =  Y^. 

q  u 

P  P  v 

Let  y^  be  given  by  an  explicit  multistep  formula  yR  =  /,  ol  yn_^  + 

q 

h  T  8.  y  .  ,  where  h  is  constant.  Assume  the  values  at  the  previous 
.L.  i  n-i 
1=1 

steps  are  known  exactly,  that  is,  yn_^  =  Y  (t^  p  and  yn_^  =  Y(tn_p, 

P  P 

i=l,...,q.  Then  (Y C^n)  -  i-s  a  constant,  independent  of  the 

values  for  both  n  and  a.  This  can  be  seen  simply  by  substitution  into 
the  above  expression. 

In  general,  we  cannot  prove  that  systematic  errors  occur.  However, 
in  practice  this  modification  does  lead  to  a  noticeable  improvement  in 
the  accuracy  of  the  predictor. 

The  key  step  in  the  algorithm,  the  determination  of  the  set  of 
indices,  L,  for  which  the  diagonal  approximation  (14)  is  valid,  depends 
on  a  user  supplied  error  tolerance,  e.  (Error  control  is  implemented 
on  a  per  step  basis.) 

Since  h  at  the  start  is  made  quite  small,  it  is  reasonable  to  solve 
for  all  the  gRi  by  (14).  Thus,  in  this  version  of  the  code,  all  indices 
1,  2,  ...,  N  are  put  into  L  at  the  beginning  of  the  integration.  (So 
that  the  entire  Jacobian  does  not  have  to  be  evaluated,  the  diagonal 
elements,  Jn^k  ynk,  are  computed  in  a  separate  subroutine.)  Next  we 
define  a  convergence  error,  E  1,  for  each  component.  If  an  absolute 
error  criterion  is  used 


^  i  |  C2i  P2i | 
E  =  y  -  y 
c  !'n  'n  1 


(28a) 


if  a  relative  error  criterion  is  used 


i  ^  ,  C2i  P2Y  C2i  P2i. » 

EC  =  2  I (yn  -  yn  )/(yn  +  yn  >1’ 


(28b) 


where  ynP2  and  ynC2  are  found  from  (17)  and  (11),  respectively.  In 
addition  an  overall  root  mean  square  (rms)  convergence  error,  E^,  is 
defined  by 
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(29) 


E_  =  (  l  (E  i) 2 /N) 1/2  . 
i=l  C 


C2 

If  Ec  <  e  we  accept  yn  as  our  final  value  for  yn.  Otherwise  the 
Newton  method  is  repeated.  Also,  if  Eci  <  e/5,  the  ith  index  remains 
in  L.  Otherwise  the  corresponding  equation  is  solved  using  the  analogue 
of  (16)  for  subsequent  steps.  The  value  e/5  is  heuristic;  but,  it 
expresses  the  basic  idea  that  very  good  accuracy  is  required  to  continue 
solving  an  equation  by  the  diagonal  approximation. 

Normally,  as  h  increases,  the  matrix  An  becomes  less  diagonally 
dominant.  Thus,  as  the  integration  proceeds,  fewer  of  the  equations  are 
solved  using  (20) .  However,  it  is  possible  for  an  equation  to  become 
diagonally  dominant  during  the  integration.  There  is  no  easy  way  to 
detect  this  event  when  it  is  solved  as  part  of  system  (16) .  But  a  crude 
check  is  made  by  monitoring  Ec .  The  condition  Ec  <  e/1000  shows  extreme 
accuracy,  so  for  the  next  step  all  the  indices  are  again  placed  into  L. 
This  allows  a  new  determination  of  the  set  of  equations  to  be  solved 
by  the  diagonal  approximation. 

The  last  important  part  of  the  integrator  is  the  algorithm  for  con¬ 
trolling  step  size.  This  is  based  on  an  estimate  of  the  truncation  error 
and  on  the  convergence  error  E^,  defined  by  (29). 

For  simplicity,  we  use  the  truncation  error  form  for  constant  step 
size,  that  is, 


E  1  =  0.0966  h4  (t)  . 


The  fourth  derivative  of  Y1  is  approximated  by  interpolating  a  fourth 
degree  polynomial  through  the  values  yn,  yn_i,  Yn-2  >  Yn>  and  Yn-1-  The 
fourth  derivative  is  then  given  by  4.'  times  the  leading  coefficient  of 
this  polynomial.  An  rms  value  for  the  fourth  derivative  is  given  by 

{l  *(4)V/n)  1/2- 

We  then  define  h^  by  the  relation 

0.0966  h  4  Y  (-4')  =  e  .  (30) 

T  rms 
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The  quantity  h,p  is  an  upper  limit  for  the  next  step  size. 

The  step  size  also  depends  on  the  convergence  error  E^.  We  define 


hc  =  hn  [1.0  +  0.1  log1()  (e/Ec)],  (31) 

and  take  the  minimum  of  the  values  h~  and  h  as  the  next  step  size  h  . . 

1  c  n+l 

Formula  (31)  provides  for  a  slow,  steady  increase  in  step  size. 

This  approach  helps  avoid  oscillations  in  the  step  size.  The  value  of 
h  can  increase  rapidly  only  if  the  predictor  and  the  corrector  are  in 
very  close  agreement. 

In  summary,  the  major  innovation  of  the  K-integrator  is  the  method 
of  solving  the  system  of  linear  equations  associated  with  the  Jacobian. 

By  using  a  diagonal  approximation,  the  effort  required  for  this  operation 
can  be  substantially  reduced.  To  increase  the  efficiency  of  this 
procedure,  more  effort  than  usual  is  invested  in  obtaining  an  accurate 
predicted  value.  A  nonorthodox  rational  predictor  is  used,  with  a 
further  heuristic  modification  based  on  the  results  obtained  in  the 
previous  steps.  The  step-size  changing  algorithm  is  based  on  both 
truncation  and  convergence  errors.  It  attempts  to  increase  the  step 
size  by  a  small  amount  each  step  rather  than  by  a  large  amount  every 
several  steps. 


III.  COMPARISONS 

The  K-integrator  has  been  compared  with  EPISODE  for  selected 
problems.  EPISODE  can  be  run  with  several  variants;  we  use  the  backward 
differentiation  formulas  (BDF),  suitable  for  stiff  problems,  with  a 
user  supplied  Jacobian. 

Some  general  observations  can  be  made  concerning  the  two  integrators 
First,  neither  one  will  be  efficient  if  a  problem  has  eigenvalues  near 
the  imaginary  axis.  This  area  is  outside  the  region  of  stability  for 
the  K-corrector  (see  page  10)  and  the  BDF  formulas.  In  a  problem  with 
eigenvalues  -  10  +  100  i9  neither  method  performed  adequately.  Such 
examples  are  not  considered  further. 

Also,  since  both  programs  are  variable  -  step  size,  they  should 
be  able  to  handle  discontinuities,  by  reducing  the  step  size  sufficiently 
This  was  one  of  the  reasons  for  developing  EPISODE.  This  aspect  is 
tested  by  the  last  three  problems. 


9W.  H.  Enright,  T.  C.  Hull,  and  B.  Lindberg,  BIT  15  (1975),  10. 
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There  are  also  important  differences.  The  K-method  was  developed 
primarily  for  equations  involving  chemical  reactions.  These  depend  on 
reaction  rates,  most  of  which  are  known  only  approximately.  Also, 
measurements  of  concentrations  of  species  cannot  be  done  to  a  very  high 
accuracy.  So  we  normally  are  satisfied  with  one  or  two  place  accuracy, 
since  any  further  accuracy  in  solving  the  equations  will  not  be 
physically  meaningful. 

Since  the  K-integrator  is  a  fixed  order  method,  we  expect  EPISODE 
to  be  more  efficient  at  stricter  error  tolerances,  where  truncation 
error  becomes  more  important.  The  treatment  of  the  Jacobian  is  also 
important.  At  a  stricter  error  tolerance,  more  steps  will  be  taken,  and 
the  Jacobian  will  not  change  as  much  per  step.  So  EPISODE  will  in 
general  update  the  Jacobian  less  frequently  and  will  become  more  efficient. 

For  a  given  error  tolerance,  the  relative  efficiency  of  the  two 
methods  is  quite  problem  dependent.  EPISODE  is  more  efficient  if  the 
Jacobian  changes  slowly.  Otherwise,  the  form  of  the  Jacobian  is  not 
important.  The  K-integrator  will  do  better  if  the  Jacobian  is  diagonally 
dominant.  Since  it  updates  the  relevant  part  of  the  matrix  each  step, 
it  is  not  bothered  by  a  rapidly  changing  Jacobian.  In  particular,  the 
K-integrator  will  tend  to  be  more  efficient  at  the  start  of  an  integration, 
when  h  is  small  and  at  least  some  of  the  Y  values  change  rapidly.  Since 
many  stiff  problems  are  integrated  until  a  steady  state  is  attained, 

EPISODE  will  be  more  efficient  at  the  end  of  the  integration,  when  the 
Y  values  change  slowly. 

If  the  strategies  of  both  the  K  -  and  Gear  integrators  could  be 
successfully  combined,  the  result  might  lead  to  a  more  efficient,  general, 
stiff  O.D.E.  algorithm. 


IV.  NUMERICAL  RESULTS 

Both  the  K-integrator  and  EPISODE  have  been  run  for  a  set  of  twelve 
problems,  given  in  Appendix  A.  The  results  of  the  comparisons  are  in 
Tables  I  through  IV.  All  the  runs  were  made  on  a  CDC  7600  in  single 
precision. 

We  are  interested  in  the  efficiency  and  accuracy  of  the  methods. 
Efficiency  can  be  measured  by  the  run  time.  The  times  given  in  the 
tables  are  the  actual  times  used  by  the  core  integrator.  The  tables 
also  give  the  number  of  derivative  calls,  the  number  of  matrix  factor¬ 
izations,  the  average  dimension  of  the  matrices  involved,  and  the 
total  number  of  integration  steps. 

For  determining  the  error,  the  "correct"  answer  Y  is  found  using  an 
error  criterion  of  10-^  or  10“10,  depending  on  the  problem.  The  error 
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Table  I.  Error  Tolerance  is  10 


RUN  FINAL  MAX 

TIME  ERROR  ERROR 


DER  MATRIX  AVERAGE 

CALLS  FAC  DIM 


E  .039  .10 

A  K  .017  .00 


.18  75  32  10.0 

.03  88  0  0.0 


E  .022  .19 

K  .022  .00 


1.21  83  20  6.0 

.98  130  39  2.0 


E  .041 

K  .048 


.07  12.19  220  28 

.09  17.93  394  152 


4.0 

1.0 


E  .003  .46 

K  .008  .00 


.46  18  8  3.0 

.00  86  10  1.0 


.002 

* 


.01 

* 


.01 


4.0 

* 


E  .088  .16 

b  K  .049  .58 


.16  251  71  7.0 

.58  263  62  2.3 


E  .086  .00 

“  K  .079  .00 


.21  128  54  9.0 

.01  203  62  5.6 


E  14.604  .12 

K  2.062  .00 


5.41 

3.13 


4323 

684 


1318 

254 


26.0 

16.9 


E  14.870  .52 

K  10.662  .14 


3.20  883  144  64.0 

3.18  832  278  30.7 


*See  text. 
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NO.  OF 
STEPS 


57 

43 


54 

56 


117 

178 


11 

42 


6 


161 

130 


84 

99 


1938 

318 


522 

412 


jtn  cfl  ^  ^  rn  ^  tn  ^cd  ^:tn 


TABLE  II.  ERROR  TOLERANCE  IS  10 


RUN 

TIME 

FINAL 

ERROR 

MAX 

ERROR 

DER 

CALLS 

MATRIX 

FAC 

AVERAGE 

DIM 

NO.  OF 
STEPS 

.082 

.41 

.41 

179 

39 

10.0 

119 

.038 

.07 

.40 

196 

0 

0.0 

97 

.049 

.35 

1.90 

184 

24 

6.0 

115 

.039 

.08 

2.35 

236 

85 

1.6 

107 

.091 

.56 

13.05 

452 

70 

4.0 

248 

.123 

.17 

65.10 

1012 

413 

1.0 

451 

.008 

.18 

.64 

48 

13 

3.0 

28 

.011 

.05 

.05 

108 

19 

1.6 

53 

.002 

.67 

.67 

8 

6 

4.0 

6 

.016 

.00 

!oo 

116 

34 

2.7 

57 

.134 

.85 

.85 

369 

94 

7.0 

225 

.094 

3.34 

3.34 

499 

145 

1.4 

247 

.152 

.00 

.27 

267 

65 

9.0 

172 

.119 

.00 

.14 

288 

97 

5.8 

142 

75.586 

.00 

8.76 

15096 

7995 

26.0 

8934 

3.913 

.05 

10.00 

1510 

504 

14.8 

726 

17.015 

2.98 

11.22 

1337 

137 

64.0 

857 

21.911 

.80 

13,72 

2022 

638 

27.8 

1004 
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TABLE  III.  ERROR  TOLERANCE  IS  10 

RUN  FINAL  MAX  DER  MATRIX  AVERAGE  NO.  OF 

TIME  ERROR  ERROR  CALLS  FAC  DIM  STEPS 


E  .156  .23 

A  K  .106  2.12 


.74  373 

2.12  550 


42  10.0  231 

0  0.0  274 


.111  .07  20.15  416 

.095  1.20  9.44  606 


35  6.0  268 

152  1.6  287 


168 

.15 

5.49 

798 

87 

4.0 

494 

337 

.68 

166.79 

2810 

1041 

1.0 

1318 

014 

.59 

1.93 

88 

15 

3.0 

52 

016 

.47 

.47 

166 

39 

1.9 

80 

.005  1.82  1.82  23 

.019  .19  .23  116 


9  4.0  14 

34  2.7  69 


.224  1.06 

.225  12.42 


1.06  613 

12.42  1314 


120  7.0  410 

183  1.5  654 


E  .260  .04 

G  K  .288  .00 


1.11  502 

.48  728 


79  9.0  310 

224  5.9  359 


E  46.046 
K  7.333 


.00  19.49  10795 

.51  29.23  3861 


4461  26.0  6533 

952  12.8  1886 


E  25.651  1.19 

K  394.610  .71 


7.00  2306  177  64.0  1688 

44.62  15719  6578  43.5  7813 
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TABLE  IV.  ERROR  TOLERANCE  IS  10 


RUN  FINAL  MAX  DER  MATRIX  AVERAGE 

TIME  ERROR  ERROR  CALLS  FAC  DIM 


19.650 

.23 

3.49 

889 

208 

64.0 

20.159 

.01 

.75 

976 

367 

37.7 

61.553 

.43 

.43 

2214 

713 

64.0 

28.872 

.02 

.23 

2322 

907 

27.1 

38.228 

.16 

8.68 

1313 

455 

64.0 

5.695 

.01 

2.87 

779 

250 

22.3 

NO.  OF 
STEPS 


525 

457 


1376 

901 


808 

361 
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may  be  absolute  (  £  (Y^y1)  /N)  ^  or  relative  (  £  ( (Y1-y:L ) /Y1)  /N)  ^  . 

i=l  i=l 

Error  is  measured  at  four  times,  equally  spaced  on  a  logarithmic  scale 
between  the  initial  step  size  hQ  and  the  final  time  tf. 

Both  methods  actually  control  error  on  a  per  step  basis.  The  error 
that  we  have  computed  is  the  global  error,  which  is  of  more  interest  to 
the  user.  The  relation  between  the  user  supplied  local  error  tolerance 
e  and  the  global  error  is  highly  problem  dependent.  However,  a  reliable 
code  should  keep  global  errors  below  a  bound  proportional  to  the  local 
error  tolerance,  for  a  given  problem.  The  tables  give  the  error  at  the 
end  time  and  the  maximum  error  at  the  four  output  times.  For  convenience, 
the  errors  are  given  in  units  of  the  error  tolerance. 

Appendix  B  gives  a  listing  of  the  code  for  the  K-integrator  for  a 
sample  problem.  The  output  is  given  at  different  times  and  for  different 
error  criterions,  as  described  above. 


The  first  five  problems  are  from  an  article  by  Enright,  Hull,  and 
Lindbergh.  These  are  part  of  a  larger  set  of  problems,  proposed  as  a 
test  set  for  comparing  stiff  integrators.  All  the  systems  are  small, 
but  they  show  that  both  methods  can  handle  a  variety  of  types  of 
problems.  All  the  problems  are  run  with  an  absolute  error  criterion. 


The  K-method  is  very  efficient  for  problems  A  and  B,  because  it  can 
take  advantage  of  the  diagonal  dominance  of  the  Jacobians.  In  fact,  in 
solving  problem  A  no  matrix  factorizations  are  performed. 


For  problem  C,  the  K-integrator  becomes  less  efficient  and  less 


-6 


The  difficulty  is  due  to  the  fact  that  one 


accurate  at  e  =  10 
component  becomes  very  large,  on  the  order  of  104.  Since  an  absolute 
error  criterion  is  being  used,  extreme  accuracy  is  required.  For  this 
type  of  problem,  the  lack  of  higher  order  formulas  in  the  K-method  can 
cause  difficulties. 


For  problems  D  and  E,  EPISODE  is  more  efficient.  The  difference 
here  is  due  to  EPISODE !s  step  changing  strategy.  Both  problems  are 
relatively  easy,  so  EPISODE  increases  the  step  size  rapidly.  The  K- 
method  has  a  more  conservative  algorithm  for  changing  step  size,  and 
takes  many  more  steps.  However,  both  programs  are  extremely  fast  here. 

_2 

The  K-integrator  does  fail  on  problem  E  for  e  =  10  .  Here  the 

problem  is  caused  by  the  fact  that  the  largest  component  is  of  the 
order  10~5.  So  a  step  can  be  accepted  by  the  integrator  even  with  no 
significant  digits  in  the  final  y  values.  As  a  result,  the  program 
uses  the  diagonal  approximation  too  long.  The  resulting  inaccuracies 
lead  to  instabilities  that  do  not  damp  out. 
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In  general,  both  methods  are  efficient  and  reasonably  accurate  on 
these  problems.  However,  the  K-integrator  does  experience  some  diffi¬ 
culties,  due  to  very  large  or  very  small  components  combined  with  an 
absolute  error  criterion,  whereas  EPISODE  does  not. 

The  remaining  systems  are  chemical  reaction  problems,  and  here  a 
relative  error  criterion  is  used.  To  prevent  control  of  the  step  size 
by  species  whose  concentrations  have  become  negligible,  an  artificial 
formation  term  of  10~30  is  added  to  y  every  time  the  derivative  sub¬ 
routine  is  called.  This  is  also  an  easy  way  of  preventing  underflow. 

Problem  F  is  a  demonstration  reaction  set  proposed  by  Edelson10. 

It  is  a  simplified  version  of  an  atmospheric  chemistry  problem.  The 
K-method  is  faster  than  EPISODE  for  e  =  10"2  and  e  =  10"4.  The  problem 
is  sufficiently  diagonally  dominant  that  the  K-integrator  can  just  use 
the  diagonal  approximation  for  a  good  part  of  the  integration,  and  works 
with  fairly  small  matrixes  even  near  the  end  of  the  integration  time. 
However,  at  e  =  10~6  its  run  time  is  equal  to  EPISODE,  and  it  is  less 
accurate.  Again,  the  K-method  will  sometimes  experience  difficulties  at 
e  =  10“  ,  due  to  the  lack  of  higher  order  formulas. 

Problem  G  and  H  are  simulations  of  the  chemistry  in  a  gun  barrel, 
i.e.,  under  conditions  of  high  temperature  and  high  pressure.  The 
problem  is  quite  stiff. 

Problem  G  involves  nine  species.  Both  programs  are  roughly  com¬ 
parable. 

Problem  H  involves  twenty-six  species,  some  of  whose  concentrations 
are  relatively  small.  EPISODE  experiences  major  difficulties,  due  to 
its  step  size  changing  strategy.  EPISODE  attempts  to  make  large  changes 
in  the  step  size,  on  the  order  of  33%.  This  leads  to  instability  in 
the  minor  species.  The  program  cannot  meet  its  error  criterion  and 
must  reduce  the  step  size.  This  prevents  the  step  size  from  increasing 
normally,  and  the  integration  takes  a  very  large  number  of  steps.  The 
K- integrator ,  with  changes  of  5%  to  10%,  can  increase  its  step  size 
consistently. 

Problem  II  is  an  atmospheric  model  of  charge  flow  under  the 
influence  of  a  large  electron  flux.  The  electron  density  starts  at  a 
high  level  and  decays  to  zero.  The  systems  consist  of  64  species. 

The  reactions  involved  are  given  by  Heimerl  and  Niles* 11. 

_2 

At  e  =  10  ,  the  K-integrator  is  somewhat  faster.  However,  the 

integrators  proceeded  in  very  different  ways.  For  example,  to  reach  the 
third  output  time,  3.16  seconds,  takes  EPISODE  12.4  seconds  and  the  K- 
integrator  only  4.5  seconds. 


10 

D.  Edelson ,  J.  Chem.  Ed.  52  (1975),  642. 
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The  K-method  is  extremely  efficient  at  the  start  of  the  integration, 
where  it  works  with  very  small  matrixes.  But  by  the  end  of  the 
integration,  the  K-integrator  is  working  with  fifty-four  by  fifty-four 
matrixes.  Meanwhile,  EPISODE  updates  the  Jacobian  about  every  3.6 
steps,  and  is  reasonably  efficient  throughout  the  integration.  So  the 
efficiency  of  the  programs  is  highly  problem  dependent. 

At  e  =  10  the  K-integrator  is  very  slow.  Because  of  truncation 
error,  it  cannot  increase  the  step  size  adequately. 

The  remaining  three  problems  involve  discontinuities  in  the  driving 
function,  that  is,  the  electron  density,  for  the  above  problem.  The 
discontinuities  occur  at  the  powers  of  10,  starting  at  10“6  seconds. 

For  comparison  purposes,  the  equations  are  integrated  out  to  a  time 
between  discontinuities,  since  the  values  are  somewhat  ambiguous  at  the 
discontinuities.  Results  are  given  only  for  e  =  10"^. 

12  is  a  sawtooth  driving  function.  The  K-integrator  is  slightly 
less  efficient  here. 

13  is  a  square  wave  driving  function.  The  comparison  is  made  at 

5  x  10^  because  EPISODE  cannot  integrate  past  the  discontinuity  at  10^. 
Even  using  the  smallest  increment  possible  on  a  CDC  7600  in  single 
precision,  the  program  cannot  meet  its  error  criterion. 

-12 

The  program  enters  an  infinite  loop.  Using  a  step  size  of  4  x  10 
it  attempts  to  cross  the  discontinuity.  Failing  its  convergence  tests 
it  reduces  the  step  size  to  2  x  10“12,  But  10^  +  2  x  10“12  ~  10^,  and 
no  progress  is  made.  It  increases  the  step  size  to  3  x  lO-!^,  which  is 
still  too  small  an  increment  to  register.  Finally,  it  increases  to  the 
value  4  x  10“^  again.  This  process  repeats  indefinitely. 

The  reason  for  this  is  the  way  EPISODE  determines  the  step  size. 
EPISODE  compares  the  original  predictor  with  the  final  accepted  y  value. 
At  a  discontinuity,  the  original  explicit  predictor  behaves  very  badly, 
since  it  has  no  information  concerning  the  changed  conditions.  The  step 
size  must  be  reduced  severely  to  obtain  agreement. 

P2 

The  K-integrator,  however,  compares  the  modified  predictor  y  with 
the  corresponding  corrector  y^2.  The  existence  of  the  discontinuity  is 
fed  in  through  the  Jacobian,  and  the  step  size  does  not  have  to  be 
reduced  so  far. 

The  K-method  is  also  much  faster  than  EPISODE.  At  the  discontinu¬ 
ities,  where  h  is  reduced,  it  can  use  the  simpler  diagonal  approximation 
to  solve  the  equations. 

Problem  14  is  a  combination  of  12  and  13.  The  electron  density 
forms  a  series  of  ramps,  alternately  increasing  and  decreasing,  but 
starting  at  the  same  value  at  each  decade.  So  there  is  a  discontinuous 
derivative  as  well  as  function.  ■ 
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Here  EPISODE  cannot  get  past  the  discontinuity  at  1.0  second;  so 
the  comparison  was  made  at  0.5  second.  Not  only  is  the  K-method  able 
to  integrate  the  problem  more  rapidly  than  EPISODE,  but  it  also  is 
able  to  integrate  over  the  entire  range  of  interest,  i.e.,  10^  seconds. 

In  conclusion  the  K-integrator  does  appear  to  be  comparable  to  the 
variable-order  EPISODE  program,  at  least  for  the  looser  error  tolerance. 
It  can  be  more  efficient,  depending  on  the  problem,  and  permits  integra¬ 
tion  over  some  severe  functional  discontinuities. 
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APPENDIX  A.  SPECIFICATION  OF  PROBLEMS 


Listed  below  are  the  systems  of  differential  equations  used  in  the 
tests.  In  all  cases  the  initial  time  tg  =  0.  The  final  time  tf  and  the 
initial  step  size  hg  are  given.  The  last  four  systems  are  not  specified 
in  detail  because  of  their  complexity. 


A. 


Linear  with  real  eigenvalues 
Y1  =  i5  Y1  Y1 (0)  =1  i  =  1,2, 


10 


t 


f 


1 


10 


-5 


B.  Linear  with  non-real  eigenvalues. 

'1  1  2 

Y  =  -10Y  +  25  Y 

'2  1  2 

Y  =  -25Y  -  10  Y 

‘3  3 

Y  =  -  4Y 

*4  4 

Y  =  -  Y 

Y5  =  -  0.5  Y5 

Y6  =  -  0.1  Y6 

tf  =20  hQ  =  10-2 


Y1  (0)  =  1 
Y2 (0)  =  1 
Y3(0)  =  1 
Y4  CO)  =  1 
Y5  CO)  -  1 
Y6  CO)  =  1 


c. 


Non- 

linear  coupling. 

Y1  - 

-  Y1 

+  2 

Y1  (0) 

=  1 

Y2  = 

-  10 

Y2  + 

20 

(Y1)2 

Y2  CO) 

=  1 

?  = 

-  40 

Y3  + 

80 

[ CY1) 2  + 

CY2)2] 

y3(0) 

=  1 

• 

II 

-100 

Y4  + 

200 [ C Y 1 ) 2  + 

CY2)2  + 

Cy3)2]  y4 

(0) 

20 

ho  =  10" 

2 
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D.  Non-linear  with  real  eigenvalues. 

*  1  1  O  7  i 

Y  =  -  Y  +  10  Y*'  (1  -  Y  ) 

'  2  2  7  ^  2 

Y  =  -  10  Y  +  3  x  10  Y  (1  -  Y  ) 

'3  '1  '2 

Y  =  -  Y  -  Y 

tf  =  1  h0  =  3-3  x  10~ 


Y1 CO)  =  1 
Y2(0)  =  0 
Y3(0)  =  0 


E.  Non-linear  with  non-real  eigenvalues. 


Y1  =  -  7.89x10 

10  yl 

7  13 

-  1.1x10  Y  Y 

Y1 (0)  =  1.76x10  3 

‘2  -10 
Y  =  7.89x10 

Y1 

9  2  3 

-  1.13x10  YY 

Y2(0)  =  6 

”3  -10 

Y  =  7.89x10 

Y1 

7  13 

-  1.1x10  Y  Y 

Y3(0)  =  0 

+  1.13xl03 

'4  7  1 

Y  =  1.1x10  Y 

Y4  - 

Y3  - 

9  2  3 

1.13x10  YY 

3  4 

1.13x10  Y 

Y4  CO)  =  0 

=  1000  hQ  =  5  x  10 


F.  Atmospheric  reaction  set.  A  set  of  7  chemical  species  and  7  reactions. 
tf  =  1000  hQ  =  10'16 


G.  Hot  Gas.  A  set  of  9  chemical  species  and  57  reactions  under  high 
temperature  and  pressure. 

tf  =  10"5  hQ  =  10-15 


H. 


Expanded  Hot  Gas.  A  set  of  26  chemical  species  and  227  reactions 
under  high  temperature  and  pressure. 


10 


-15 
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I.  BENCHMARK  -76.  A  set  of  64  chemical  species  and  498  reactions 
simulating  chemistry  in  the  upper  atmosphere.  The  reactions  are 
driven  by  the  electron  density. 


Relaxation  of  driving  force. 

,  =  104  -10 
*f  ho  ■  10 


2.  Saw  tooth  drawing  force  (discontinuous  derivative). 


t^.  =  5  x  10" 


h0  =  10 


-10 


3.  Square  Wave  driving  force  (discontinuous  function) . 


t^  =  5  x  10 


h0  =  10 


-10 


4.  Ramp  driving  force  (discontinuous  function  and  derivative) 

-10 


t^.  =  0.5 


h0  -  10 
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APPENDIX  B.  PROGRAM  LISTING 

A  listing  of  the  computer  code  follows.  The  code  is  set-up  to 
solve  problem  C,  page  27. 
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Please  take  a  few  minutes  to  answer  the  questions  below;  tear  out 
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Maryland  21005.  Your  comments  will  provide  us  with  information 
for  improving  future  reports. 

1 .  BRL  Report  Number _ 

2.  Does  this  report  satisfy  a  need?  (Comment  on  purpose,  related 
project,  or  other  area  of  interest  for  which  report  will  be  used.) 


3.  How,  specifically,  is  the  report  being  used?  (Information 
source,  design  data  or  procedure,  management  procedure,  source  of 
ideas,  etc.) _ 


4.  Has  the  information  in  this  report  led  to  any  quantitative 
savings  as  far  as  man-hours/contract  dollars  saved,  operating  costs 
avoided,  efficiencies  achieved,  etc.?  If  so,  please  elaborate. 


5,  General  Comments  (Indicate  what  you  think  should  be  changed  to 
make  this  report  and  future  reports  of  this  type  more  responsive 
to  your  needs,  more  usable,  improve  readability,  etc.) 


6.  If  you  would  like  to  be  contacted  by  the  personnel  who  prepared 
this  report  to  raise  specific  questions  or  discuss  the  topic, 
please  fill  in  the  following  information. 


Name: 
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